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ABSTRACT 

We use a suite of numerical simulations to investigate the mechanisms and effects of radial 
migration of stars in disk galaxies like the Milky Way (MW). An isolated, collisionless stellar 
disk with a MW-like scale-height shows only the radial "blurring" expected from epicyclic or- 
bits. Reducing the disk thickness or adding gas to the disk substantially increases the level of 
radial migration, induced by interaction with transient spiral arms and/or a central bar. We also 
examine collisionless disks subjected to gravitational perturbations from a cosmologically 
motivated satellite accretion history. In the perturbed disk that best reproduces the observed 
properties of the MW, 20% of stars that end up in the solar annulus 7 kpc < R < 9 kpc 
started at R < 6 kpc, and 7% started at R > 10 kpc. This level of migration would add 
considerable dispersion to the age-metallicity relation of solar neighborhood stars. In the iso- 
lated disk models, the probability of migration traces the disk's radial mass profile, but in 
perturbed disks migration occurs preferentially at large radii, where the disk is more weakly 
bound. The orbital dynamics of migrating particles are also different in isolated and perturbed 
disks: satellite perturbations drive particles to lower angular momentum for a given change 
in radius. Thus, satellite perturbations appear to be a distinct mechanism for inducing radial 
migration, which can operate in concert with migration induced by bars and spiral structure. 
We investigate correlations between changes in radius and changes in orbital circularity or 
vertical energy, identifying signatures that might be used to test models and distinguish radial 
migration mechanisms in chemo-dynamical surveys of the MW disk. 

Key words: methods: numerical; Galaxy: kinematics and dynamics; Galaxy: disc; galaxies: 
formation; galaxies:evolution 



1 INTRODUCTION 

Under the influence of gravity, disk galaxies are expected to assem- 
ble in an "inside-out" fashion: stars form first from high-density gas 
in the central region of the galaxy where the potential is deepest, 
and s ubsequently at increasing galacto-centric radii (e.g. iLarsonl 
fl976h . An immediate consequence of this formation scenario is 
that stars bom at the same time and in the same region of a galaxy 
should have similar chemical compositions. However, observations 
in our Galaxy suggest that these initial conditions are not main- 
tained. IWielen et al.l jl996l) argued that the Sun was substantially 
more metal rich than nearby solar age stars and the local interstellar 
medium (ISM). A recent recalibration of the Geneva Copenhagen 
Survey (GCS) using th e infrared flux method finds no discrep- 
ancy with solar age stars dCasagrande et al.l201 ll) . but even modern 
studies confirm that the age-metallicity relationships (AMRs) of 
field and solar neighborh ood stars are characterized by higher dis- 
persions than expecte d l Edvardsson et alj 1993 : Nordstrom et al.l 
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simple chemical evolution models that divide the Galaxy into in- 
dependently evolving concentric annuli predict many more low 
metallicity G-dwarfs in our region of the disk compared to those 
observed, a discrepancy known as "the local G-dwarf problem" 
dvan den Berghll 19621 : ISchmiddfl963h . Evidence seemingly in con- 
tradiction to standard galaxy chemical evolution theory is not lim- 
ited to our own Galaxy. Metallicity gradients in dis k galaxies are 
shallower than predicted by class ical models (e.g., Magrini et al] 
120071) . iFerguson & Johnson! d200ll) and lFerguson et al.l d2007h find 
unexpectedly old stellar populations on nearly circular orbits in the 
outskirts of M31 and M33, respectively. The outermost regions of 
NGC300 and NG C7739 show flattened or p ositive abundance gra- 
dients with radius dVlaiic etal . 20091 1201 Jj). These perplexing ob- 
servations cannot be readily explained within the confines of classic 
galaxy formation models. 

A natural explanation for the observational challenges above 
arises if the present day radii of many stars could be signifi- 
cantly different from their birth radii. One difficulty in establish- 
ing radial migration as a common phenomenon, from a dynamical 
standpoint, lies in finding a mechanism that can cause a substan- 
tial fraction of stars to migrate several kiloparsecs while retain- 
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ing the observed approxim ately circular orbits. In a seminal paper, 
ISellwood & Binnevl J2002L hereafter SB02) investigated the rela- 
tionship between changes in stellar angular momentum and disk 
heating. They found that radial migration is a ubiquitous process in 
spiral galaxies; stars naturally migrate (change angular momentum) 
as they resonantly interact with transient spiral waves. Stars in coro- 
tational resonance (CR) with said waves are scattered without heat- 
ing the disk and maintain their nearly circular orbits (unlike Lind- 
blad resonance (LR) scattering). In the present paper, we examine 
radial migration in simulations of disk galaxies. Our experiments 
include galactic disks evolved both in isolation and under the action 
of infalling satellites of the type expected in the currently favored 
cold dark m atter (CDM) paradigm of hierarchica l structure for- 
mation (e.g.. IPeeblesIl 19821 : iBlumenthal et alJI 19841). The latter set 
of exp eriments were prese nted in the studies of Kazantzidis et al.1 
J2008L hereafter K08) and lKazantzidis et alj |2009|) and were uti- 
lized to investigate the generic dynamical and morphological signa- 
tures of galactic disks subject to bombardment by CDM substruc- 
ture. 

Inspired by SB02, several groups have recently investigated 
the potential role of radial mixing in the chemical and dynamical 
evolution of disk galaxies. ISchonrich & Binnevl d2009l) presented 
the first chemical evolution model to incorporate radial migration. 
The rate at which stars migrate via the SB02 mechanism is left as a 
free parameter constrained by the metallicity dist ribution function 
(MDF ) of solar neighborhood stars in the GCS dNordstrom et al.1 
120041) . Their model successfully reproduced, within systematic un- 
certai nties, the observed age -metallicity distribution of stars in the 
GCS dHolmberg et"ai]|2007l) and the observed correlation between 
tangen tial velocity and abundance pattern described by iHavwoocj 
J2008I) . However, there is partial degeneracy between the magni- 
tude of radial migration and other parameters in the model such as 
star-formation rates and gas inflow characteristics. Furthermore, it 
is unclear whether the level of migration required to fit the data is 
consistent with theoretical expectations . 

Numerical simulations have confirme d the occurrence of r a- 
dial migration under a variety of conditions. iRoskar et al.l d2008al lbT) 
studied the migration of stars in a simulation of an isolated 
Milky Way (MW)-sized stellar disk formed from the cooling of 
a pressure-supported gas cloud in a 1O 12 M0 dark matter halo. 
In their simulations, some older stars radially migrated to the 
outskirts of the disk while maintaining nearly circular orbits, 
forming a population akin to that observed in M31 and M33 
( [Ferguson & Johnson] l200ll : iFerguson et al.l |2007|) . IRoskar et al.1 
d2008bl) found that ~ 50% of all stars in the solar neighborhood 
were not born in situ; this is a natural explanation for the observed 
dispersion in the AMR and solar neighborhood metallicity distri- 
bution function (MDF). 

More recently. lOuillen et al.ld2009T) investigated radial migra- 
tion in a stellar disk perturbed by a low-mass (~5x 10 9 M ) 
orbiting satellite. Their numerical simulations integrated test par- 
ticle orbits in a static galactic potential and highlighted the fact 
that mergers and perturbations from satellite galaxies and subha- 
los can induce stellar radial mixing. Although informative, test par- 
ticle simulations in a static isothermal potential will not capture 
all the relevant physics of the process of stellar radial migration 
in disk galaxies; the interactions between the gravitational pertur- 
bations and the self-gravity of the disk are essential to a detailed 
analys is of the phenomeno n. In our paper, we expand on the anal- 
ysis of lOuillen et"all d2009l) by investigating radial migration using 
fully self-consistent numerical simulations both with and without 
satellite bombardment. 



There are now several established phenomena that can cause 
a star to populate a region of the disk different from its birth radii. 
Stars on elliptical orbits maintain their guiding center and angu- 
lar momentum (modulo asymmetries in the potential) but can be 
found over the range in galacto-centric radius defined by their peri- 
center and apocenter. Changing a star's angular momentum, and 
hence its guiding center, requires direct scattering or a resonant 
interaction with transient patterns in the disk. The local encoun- 
ters o f stars with molecular clouds (e.g. ISpitzer & Schwarzschild 
Il953l) or Lindblad resonance (LR) scattering between stars and spi- 
ral waves (e.g. SB02) both change stellar guiding centers (albeit 
to a relatively small degree) and increase the random motions of 
stars over time, "blurring" the disk. Stars scattered at CR with spi- 
ral waves can change their guiding centers by several kiloparsecs 
without increasing the amplitude of their radial motion. For any 
single spiral wave, SB02 predict that stars are scattered on each 
side of the CR, "churning" the contents of the dislQ. Stars may 
undergo several encounters with transient spiral waves throughout 
their lifetimes. While SB02 investigate all resonant interactions be- 
tween spiral waves and stars, we will refer to this special case of 
CR as the "SB02 mechanism'0. As SB02 note, migration due to 
spiral waves can be described by blurring and churning regardless 
of how the waves arise (satellites, e.g., could induce spiral struc- 
ture that would lead to migration described by SB02). Simulations 
have shown that other transient wave patterns internal to a galaxy, 
such as those resulting from bar propagation, can produce reso- 
nance o verlap with existing spiral patterns and induce radial mi- 
gration dMinchev & Famaey||201ol : iMinchev et alll201ll) . Orbiting 
satellites, external to the galaxy and discussed above, will have a 
complex interaction with the disk as they provide a means of di- 
rect scattering over a large area and also induce spiral modes in the 
disk. In this work, we aim to characterize radial migration induced 
by satellite bombardment and compare its effects on the stellar disk 
with those observed in secularly evolved galaxies. 

Our investigation complements earlier and ongoing radial mi- 
gration studies. We perform a simulation campaign, including nu- 
merical experiments of isolated disk galaxies with different scale 
heights and gas fractions, which in turn lead to different levels of 
spiral structure. For the first time, we examine the effect of satel- 
lite bombardment on radial migration utilizing simulations where 
galactic disks are subjected to a cosmologically motivated satellite 
accretion history. Via a comparative approach, we determine how 
the magnitude and efficiency of radial migration depend on input 
physics, establish correlations between orbital parameters and mi- 
gration, and present evidence that each of the three migration mech- 
anisms is distinct in the examined parameter space. These charac- 
teristics lead to possible observational signatures that may constrain 
the relative importance of each migration mechanism in the Milky 
Way. 



2 METHODS 

2.1 Isolated Disk Models 

We employ the method of IWidrow & Dubinskil d2005l) to con- 
struct numerical realizations of self-consistent, multicomponent 



1 Blu rring and churning are the terms proposed by ISchonrich & Binnevl 
120091) to describe these distinct aspects of radial migration. 

2 In the recent radial migration studies of secularly evolved simulations, 
significant migration is almost always attributed to CR scattering. 
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disk galaxies. These galaxy models consist o f an ex ponential stel- 
lar dis k, a Hemquist bulge dHemquisdl 19901 and a lNavarro et al.l 
d 19961 hereafter NFW) dark matter halo. They are characterized by 
15 free parameters that may be tuned to fit a wide range of obser- 
vational data for actual gal axies including the MW and M31. The 
IWidrow & Dubinskil l l2005l) models are derived from three-integral, 
composite distribution functions and thus represent self-consistent 
equilibrium solutions to the coupled Poisson and collisionless 
Boltzmann equations. Owing to their self-consistency, these galaxy 
models are ideally suited for investigating th e complex dynamics 
involv ed in the process of radial mixing. The IWidrow & Dubinskil 
< l2005t) method has been recently used in a variety of numeri- 
cal studies associated with instabilities in disk galaxies, includ- 
ing the dynamics of warp s and bars jDubinski & Chakrabartvl 
1 20091 : iDubinski et alj[2009l) . the gravit ational interaction be tween 
galactic disks and infalling satellites jGauthier e t al. 2006; K08; 
IPurcell et al.l2009l;lKazantzidis et alj2009l) . and the transformation 
of disky dwarfs to dwarf spheroi dal galaxies under the a ction of 
tidal forces fro m a massive host (Kazantzidis et"al]|201ll) . We re- 
fer the reader to lWidrow & Dubinskil d2005l) for an overview of all 
relevant parameters and a detailed description of this technique. 

For the majority of numerical experiments in the presen t 
study, we employ model "MWb" in IWidrow & Dubinskil Hoblh . 
which satisfies a broad range of observational constraints on the 
MW galaxy. Specifically, the stellar disk has a mass of A/disk = 
3.53 x 10 10 M Q , a radial scale length of R d = 2.82 kpc, and 
a sech 2 scale height of Zd = 400 pc. We note that the adopted 
value for the scale height is consiste nt with that inferred for the 
old, t hin stellar disk of the MW (e.g-. lKent et alJI 199ll : ljuric et all 
120081) . The equivalent exponential scale height is approximately 
200 pc, but the sech 2 vertical distribution is more accurate. The 
bulge has a mass and a scale radius of Alt = 1.18 x 10 10 Mq 
and a,b — 0.88 kpc, respectively. The NFW dark matter halo has 
a tidal radius of Rh = 244.5 kpc to keep the total mass finite 
(K08), a mass of M h = 7.35 x 10 11 M , and a scale radius of 
rh = 8.82 kpc. The total circular velocity of the galaxy model 
at the solar radius, Rq ~ 8 kpc, is V c (Rq) = 234.1 kms -1 , 
and the Toomre disk stability parameter is equal to Q — 2.2 at 
R — 2.5Rd- We note that direct numerical simulations of the evo- 
lution of model MWb in isolation confirm its stability against bar 
formation for 10 Gyr. 

We wish to address the dependence of radial mixing in iso- 
lated disk galaxies upon initial disk thickness and the presence of 
gas in the galactic disk. Bec ause of their smaller "effe ctive" Toomre 
Q stability parameter (e.g.,|Romeo & WiegertEoi ]]), thinner disks 
are more prone to gravitational instabilities and thus yield stronger 
and better defined spiral structure. Therefore, we might expect them 
to cause enhanced radial mixing compared to their thicker coun- 
terparts. Similarly, because radiative cooling in the gas damps its 
random velocities, the presenc e of gas is associated with lo ng-lived 
spiral structure in disks (e.g., ICarlberg & Freedmanlll985l) . which 
may act to increase the amount of radial mixing. Correspondingly, 
we initialize several additional disk galaxy models. 

The first modified galaxy model was constructed with the 
same parameter set as MWb but with a scale height 2 times smaller 
(zd = 200 pc). Except for disk thickness and vertical velocity dis- 
persion, all of the other gross properties of the three galactic com- 
ponents of this model are within a few percent of the correspond- 
ing ones for MWb. Not surprisingly, being more prone to gravi- 
tational instabilities, this model develops a bar inside ~ 5 kpc at 
time t ~ 1.2 Gyr. The second set of modified galaxy models are the 
same as MWb except for the fact that a fraction f g of the mass of 



the initial stellar disk is replaced by gas. Thus, the resulting gaseous 
component is constructed with the same initial density distribution 
as the stellar disk. In the present study, we employ two values for 
the gas fraction, f g = 0.2 and f g — 0.4, and the adopted methodol- 
ogy will be described in detail in a forthcoming paper (Kazantzidis 
et al. 201 1, in preparation). Briefly, the construction of the gas disk 
is done by assuming that its vertical structure is governed by hydro- 
static equilibrium and by computing the potential and the resulting 
force field for the radially varying density structure of the gas com- 
ponent. By specifying the polytropic index and the mean molecular 
weight of the gas and using a tree structure for the potential cal- 
culation, the gas azimuthal streaming velocity is determined by the 
balance between gravity and centrifugal and pressure support. 



2.2 Perturbed Disk Models 

In the standard CDM paradigm of hierarchical structure formation, 
galaxies grow via continuous accretion of smaller satellite systems. 
To investigate the effect of accretion events on radial mixing in 
galactic disks, we have also analyzed iV-body simulations of the 
gravitational interaction between a population of dark matter sub- 
halos and disk galaxies. These simulations were first presented by 
K08, and we summarize them briefly here. 

K08 performed high-resolution, collisionless iV-body simula- 
tions to study the response of the galactic model MWb discussed 
above to a typical ACDM-motivated satellite accretion history. The 
specific merger history was derived from cosmological simulations 
of the formation of Galaxy-sized CDM halos and involved six dark 
matter satellites with masses, orbital pericenters, and tidal radii 
of 7.4 x 10 9 < M sat /M < 2 x 10 10 , r pcri < 20 kpc, and 
rtid > 20 kpc, respectively, crossing the disk in the past ~ 8 Gyr. 
As a comparison, the Large Magellanic Clo ud has an estimated to- 
tal present mass of ~ 2 x 10 10 M Q (e.g JSchommer et al.|[l992l; 
iMastropietro et al.l2005l) . so the masses of the perturbing dark mat- 
ter satellites correspond to the upper limit of the mass function of 
observed satellites in the Local Group. K08 modeled satellite im- 
pacts as a sequence of encounters. Specifically, starting with the 
first subhalo, they included subsequent systems at the epoch when 
they were recorded in the cosmological simulation. We note that 
although K08 followed the accretion histories of host halos since 
z ~ 1, when time intervals between subhalo passages were larger 
than the timescale needed for the disk to relax after the previous 
interaction, the next satellite was introduced immediately after the 
disk had settled from the previous encounter. Each satellite was re- 
moved from the simulation once it reached its maximum distance 
from the disk after crossing. This approach was dictated by the de- 
sire to minimize computational time and resulted in a total simu- 
lation time of ~ 2.5 Gyr instead of ~ 8 Gyr that would formally 
correspond to z = 1. 

K08 and Kazantzidis et al ] (HI) found that these accretion 
events severely perturbed the galactic disk of model MWb without 
destroying it and produced a wealth of distinctive morphological 
and dynamical signatures on its structure and kinematics. In this 
paper, we will investigate the magnitude of radial mixing induced 
in disk model MWb by these accretion events. As in the case of iso- 
lated disk models, it is worthwhile to examine the dependence of 
radial mixing upon disk thickness. For this purpose, we employed 
the thinner galaxy model with a scale height of Zd — 200 pc de- 
scribed above and repeated the satellite-disk encounter simulations 
of K08. 
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Figure 1. Surface density maps of the stellar distributions of galactic disks with different initial scale-heights z d . Maps are presented only for the collisionless 
experiments and include the initial (left panels), final isolated (middle panels), and final perturbed disks (right panels). Each panel includes the face-on (bottom 
panels) and edge-on (upper panels) distributions of disk stars. Particles are color-coded on a logarithmic scale, with hues ranging from blue to white indicating 
increasing stellar density. Local density is calculated using an SPH smoothing kernel of 32 neighbors. The density ranges from 7 X 10 5 to 2 X 10 9 Mq / kpc 2 . 



2.3 Numerical Parameters 

All collisionless numerical simulations discussed in this paper were 
carried out with the mul ti-stepping, parallel, tree iV-body code 
PKDGRAV tStade<l200ih . The hydrodynamical simulations were 
performed with the p arallel TreeSPH A-body code GASOLINE 
dWadslev et al. I l2004j) . In the gasdynamical experiments, we in- 
clude atomic cooling for a primordial mixture of hydrogen and 
helium, star formation and (thermal) feedb ack from supernovae . 
Our star formation recip e follows th at of IStinson et all fcOOfj) . 
which is based on that of lKatzl <1992h . Gas particles in cold and 
dense regions which are simultaneously parts of converging flows 
spawn star particles with a given efficiency c* at a rate propor- 
tional to the local dynamical time. Feedback fro m supernovae 
is trea ted using the blast-wave model described in IStinson et al .1 
( l2006h . w hich is based on the analy tic treatment of blastwaves de- 
scribed in lMcKee & Ostrikerl i ll 9771) . In our particular applications, 
gas particles are eligible to form stars if their density exceeds 0. 1 
atoms/cm 3 and their temperature drops below T max = 1.5 x 10 4 K, 
and the energy deposited by each Type-II supernova into the sur- 
rounding gas is 4 x 10 erg. We note that this choice of parameters 



and numerical techniques is sho wn to produce realistic disk galax- 
ies in cosmological simulations dGovernato et al ] l2007h . Lastly, in 
an attempt to investigate the effect of star formation efficiency on 
radial mixing, for each of the gas fractions f g above, we adopted 
two different values for c*, namely c* = 0.05 (which reproduces 
the slope and normalization of the observed Schmidt law in isolated 
disk galaxies) and a lower value of c* = 0.01. 

All A-body realizations of the disk galaxy models contain 
N d = 10 6 particles of m d = 3.53 x 1O 4 M in the disk; N b 
5 x 10 5 , m b = 2.36 x 10 4 M Q in the bulge; and N h = 2 x 10 6 , 
rrid = 3.675 x 10 in the dark matter halo. The gravitational 
softening lengths for the three components were set to e d — 50 pc, 
ej = 50 pc, and = 100 pc, respectively. These are "equiva- 
lent Plummer" softenings; the force softening is a cubic spline. In 
the hydrodynamical simulations of isolated galaxies, gaseous disks 
were represented with N g — 2 x 10 5 and N g = 4 x 10 s for gas 
fractions of f g = 0.2 and f g — 0.4, respectively. In these cases, 
the gravitational softening length for the gas particles was set to 
e g = 50 pc. All numerical simulations of isolated and perturbed 
disks were analyzed at 2.5 Gyr. This time-scale corresponds to ap- 
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proximately 17.5 orbital times at the disk half-mass radius. For 
evaluating the impact of accretion events, it is fairest to compare 
isolated and perturbed disks after the same amount of integration, 
but we may underestimate the total amount of mixing in the iso- 
lated (and, to a lesser extent the perturbed) models. In the future, 
we plan to evolve selected simulations for the full 8 Gyr interval 
since 2 = 1, but such simulations will require more than three 
times the computational resources used here. 

Exchange of angular and linear momentum between the in- 
falling satellites and the disk tilt the disk plane and cause the disk 
center of mass to drif t from its initial position at the origin of the co- 
ordinate frame (K08: lKazantzidis et alj2009l ). Therefore, we calcu- 
late Ar in the satellite-disk encounter experiments after removing 
the global displacement and tilt of the disk galaxy by determining 
the principal axes of the total disk inertia tensor and rotating our 
original coordinate system such that it is aligned with this tensor. 



3 RESULTS 

Figure [TJ presents surface density maps of the stellar distribu- 
tions of our four collisionless galactic disks with different initial 
scale-heights Zd- Each simulation exhibits a distinctive combina- 
tion of morphological features that could affect radial migration. 
The initial smooth disks are unstable to spiral instabilities aris- 
ing from the swing amplification of particle shot noise, an effect 
that leads to emerging spiral stru cture as seen in the final isolated 
disks (e.g. ljulian & Toomreil966l) . Owing to its smaller "effec tive" 
Toomre Q stability parameter (e.g., iRomeo & Wieger3l201 ll) . the 
Zd = 200 pc disk (hereafter, disks with this initial scale height will 
be referred to as "thin") develops prominent spiral structure and a 
strong bar. Conversely, the Zd = 400 pc disk (hereafter, disks with 
Zd = 400 pc will be referred to as "thick") has relatively little spi- 
ral structure and does not form a bar in isolation. We emphasize 
that the "thick", Zd = 400 pc disk is the one in best agreement 
with the old, thin disk of the MW, while Zd = 200 pc is too thin 
(see Section [2~T| >. 

The radial and vertical morphology of the perturbed disks 
is distinct from their isolated counterparts. Both perturbed disks 
develop bars and are characterized by prominent flaring and 
much larger scale heig hts than those of the isolated disks (K08; 
iKazantzidis et alj|2009l ). There is evidence that some spiral struc- 
ture evident in the isolated disks has been washed out in the simula- 
tions with substructure bombardment: within 10 kpc of the galac- 
tic center, local enhancements of the stellar surface density evident 
in the isolated disks are substantially more diffuse after the action 
of the infalling satellites. Throughout this section, we will inves- 
tigate how the growth and dissipation of spiral structure affect ra- 
dial migration. While we have investigated similar maps in the four 
hydrodynamical simulations, we do not show them here as they 
are qualitatively similar to their collisionless counterparts. As ex- 
pected, however, the magnitude of spiral structure increases as the 
gas fraction rises. We note that the differential effect of gas on the 
strength of spiral structure we find here is smaller than that reported 
in previous numerical inve stigations of isolated disk galaxies (e.g., 
lBarnes&Hernquisj|l996l) . This is mainly due to the effect of stel- 
lar feedback, which causes the IS M in our simulation s to become 
turbulent and multi-phase (see also lStinson et al . 2006). 



3.1 Radial Migration 

We first investigate where particles move throughout each simula- 
tion. Figure [2] shows the distribution of Ar= ri — r\ in the six 
isolated disk simulations, where n and n refer respectively to each 
particle's final and initial projected distance from the galactic cen- 
ter. Each panel contains the Ar distribution, the median |Ar| (de- 
noted Ar me d), and 80" 1 percentile | Ar | (denoted Argo) for two of 
the six isolated galaxies in our simulation suite. 

The isolated, collisionless simulations (left panel) clearly 
demonstrate that disk scale height and gas content affect the ra- 
dial migration process. The thin disk's Ar distribution (black line) 
is considerably broader than the thick disk's (gray hatch). Both 
Ar me d= 0.91 kpc and Argo— 2.03 kpc of particles in the thin 
disk are approximately double those found in the thick disk. 

The remaining panels of Figure [2] compare the Ar distribu- 
tions of the hydrodynamical simulations of our sample. We exam- 
ine four isolated disks with two initial gas fractions: f 3 = 20% 
(middle panel) and f g = 40% (right panel). The two histograms in 
each panel represent star formation efficiencies (c*) of 1% (solid 
line) or 5% (gray, hatched histogram). Gas has a strong impact on 
particle radial migration: Ar me d increases from 0.47 kpc in the 
collisionless case (all hydrodynamical simulations are of "thick" 
disks) to ~ 0.70 kpc when f g = 20% and ~ 0.80 kpc when 
f g — 40%. The extent of radial migration is less dependent on 
star formation efficiency; there are only minor differences amongst 
each pair of histograms in the middle and right panels. However, in 
both panels, Ar me d and Arso are highest when c*= 1%. In these 
four hydrodynamical simulations, the time-averaged gas content of 
the disk is directly correlated with the percentage of particles that 
migrate and their typical displacement relative to their formation 
radius. 

Decreasing the scale height of the stellar disk, increasing the 
fraction of the initial disk mass in gas, and lowering the star- 
formation efficiency all increase the midplane density of the disk. 
This enhanced density increases the coherence and self-gravity of 
spiral structure, supporting it against the dissipative natu re of in- 
dividu al particles' random motions (as described by, e.g.. lToomrel 
1 1977b . The correlation of spiral structure strength with the fraction 
of particles that migrate away from their birth radii and the median 
migratory distance traversed suggests that the SB02 mechanism, 
intimately linked with spiral waves, plays a major role in the migra- 
tion of particles in the isolated systems. The symmetric shapes (to 
within 0.4% about kpc) of the Ar distributions are also consis- 
tent with migration via the SB02 mechanism. Resonances between 
the bar and spiral structure could also influence the timescal e for 
migration in the thin disk models jMinchev & Famaey||2010l) . but 
the Zd = 400 pc models do not develop bars, regardless of whether 
they include gas. 

Figure [3] compares the Ar distributions of the isolated colli- 
sionless disks to those of their perturbed counterparts. The Ar me d 
and Argo of each distribution are labeled and color-coded to 
match the corresponding histogram. In each panel, the most sharply 
peaked histogram (red) shows the expected Ar distribution from 
epicyclic motion alone, which we compute given the particle's ini- 
tial phase space coordinates. We model the epicycle as a simple 
harmonic oscillator about the particle's guiding center radius (R g ) 
with an amplitude set by its initial radial energy (E r ) (Chapter 3, 
iBinnev & Tremainell2008h . The energy associated with the circular 
component of the particle's orbit (i? c i rc ) is set by R g , which is the 
radius of a circular orbit with angular momentum equal to the mid- 
plane component of the particle's angular momentum. The radial 
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Figure 2. The distribution of radial shift (Ar= n — n) for all disk particles in each of the six isolated disk simulations. Ar me< j and Argo specify the 
median and 80th percentile distance traveled, respectively. Gray values refer to the hatched histograms; black values are associated with the solid black line 
in each panel. We report the fraction of all particles in non- overlapping 500 pc bins of Ar. The two collisionless simulations (left panel) illustrate the effect 
of disk structure on radial mixing, i.e., there is more radial migration in the smaller scale height (2<j = 200 pc, black line) disk than in its thicker counterpart 
(z d = 400 pc, hatched histogram). The four hydrodynamical simulations are in the middle (f g = 20%) and right (/ 9 = 40%) panels. In both of these 
panels, the gray, hatched histograms represent the simulations with higher star formation efficiency (c* = 0.05); those with c*= 0.01 are shown in black. 
Both increased gas fractions and smaller star formation efficiencies yield greater radial mixing. 



energy of the orbit is E T = E to t — E 2 — Scire where E z is the ver- 
tical energy component and E to t is the total energy . As E T sets the 
amplitude of the epicycle, we can solve for the position of the par- 
ticle as function of time with respect to its guiding center radius. 
We choose two random phases of the epicycle oscillation, desig- 
nating the first as n and the second as ry, and compute Ar. The 
resulting "baseline" distribution is shown in red and is a result of 
observing the two random phases (akin to our "initial" and "final" 
snapshots) of the initially elliptical orbits. It does not take into ac- 
count the potential heating of these initial orbits and the subsequent 
increase in amplitude of their radial motion (blurring). The shape 
of the distribution does not change if we change our random seed or 
average a larger number of phases to obtain rt and n. For the iso- 
lated thick disk (right panel), the Ar distribution from full dynami- 
cal evolution is only marginally broader than the expected baseline 
distribution, suggesting that these stars are not heated with time and 
continue to follow their initial orbits. However, satellite bombard- 
ment broadens the Ar distribution dramatically, nearly doubling 
both Ar me d and Argo, necessitating guiding center modification 
of a substantial fraction of the orbits. 

For the thin disk, isolated evolution produces a much broader 
Ar distribution than the baseline distribution, demonstrating the 
impact of the bai0 and spiral structure in this more unstable disk. In 
this case, satellite bombardment only slightly increases the width of 
the Ar distribution, despite the strong impact on disk structure that 
is visually evident in FigureQ] We suspect that the small net differ- 
ence between these two histograms reflects a cancellation between 
two competing effects of satellite bombardment. Accretion events 



3 Bars can increase the eccentricity of particles, contributing to the "blur- 
ring" of the disk, and thei r potential resonance overl ap with spiral structure 
can induce migration (see Minchev & Famaey 201ol). 



heat the stellar disk and thereby suppress the development of spi- 
ral structure, thus reducing the level of SB02 migration. However, 
the accretion events also induce radial mixing directly via their dy- 
namical perturbations. The Ar histograms of both perturbed disks, 
while still approximately symmetric, are noticeably more asym- 
metric than those of the isolated disks, with 2% (3%) more particles 
moving outwards than inwards in the thin (thick) cases, compared 
to ^ 0.4% for the isolated models. 

Figure |4]presents clear evidence that satellite-induced migra- 
tion is distinct from that in the isolated experiments, acting in dif- 
ferent environments from either epicyclic blurring or spiral-induced 
churning. In each panel, solid, dashed, and dotted lines show the 
fraction of particles at each birth radius i?f or m that migrate by more 
than | Ar | = 1.0, 2.0, or 3.0 kpc, respectively. Shaded regions 
show the scaled radial surface density profile of the initial disk. For 
both isolated disks, the migration probability decreases outwards 
and approximately traces the surface density profile. This behav- 
ior is similar to that assum ed in the chemical evolution models of 
ISchonrich & Binnevl ( T2009|) . who parametrized the probability that 
a star migrates as proportional to the mass surrounding it. How- 
ever, for both perturbed disks the probability of migration is flat or 
increasing outwards beyond iZfbrm= 10 kpc, where the disk poten- 
tial weakens, and it definitely does not trace the mass distribution. 
To better understand this difference between satellite- and spiral- 
induced radial mixing, we now investigate how changes in angular 
momentum and energy are correlated with change in radius. 

3.2 Orbital Characteristics 

Total energy and angular mo mentum are the classic two- 
dimensional integrals of motion l lBinnev & Tremainell2008|) . Fig- 
ure [5] shows Lindblad diagrams for the two initial and four final 
states of the collisionless simulations, plotting particle specific an- 
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Figure 3. The Ar distribution for all disk particles in the four collision- 
less simulations. Each histogram show the fraction of particles in non- 
overlapping 500 pc bins of Ar. The isolated (hatched, gray histogram) and 
perturbed (black line) = 200 pc (left panel) and 2S<j = 400 pc (right 
panel) disks are shown. The perturbed disks show greater dispersion in Ar 
compared to isolated galaxies with the same geometry. The red histograms 
denote the expected Ar distribution from epicyclic motion alone, given the 
initial orbital configurations of the particles in each disk. Ar me< j and Argo 
are labeled and color-coded for the three distributions in each panel. 



gular momentum projected along the axis of symmetry ( J z ) versus 
specific binding energy (E). Particles are grouped into 50 distinct 
linear bins of E. We calculate the median J z in each bin. Red, yel- 
low, and orange regions connect the 68 th , 95 th , and 99 th percentile 
J z intervals (centered on each median J z ), respectively, across all 
energy bins. The 1% most discrepant particles in J, for a given E 
are plotted as individual points. Using each disk's rotation curve, 
we plot J z and E for circular orbits in the midplane of the disk 
(green line). Here, J z — v c (r) x r , where v c (r) is the circular 
velocity at radius r and E = \v'i(r) + $(r), where $(r) is the 
potential in the disk midplane at radius r. Particles on a circular 
orbit have the maximum J z allowed given their energy. 

By construction, particles are initially (left column of Fig- 
ure |5j on nearly circular orbits with a small radial velocity com- 
ponent. The red region, falling close to the circular orbit curve in 
both initial disks, confirms that <j Vr is small. Particles significantly 
displaced from the green curve in the initial states either have ex- 
treme v T or are on highly inclined orbits ( J z is projected along the z 
axis). Due to this inclination effect, the initial thick disk has a lower 
median J z as a function of E than the thin disk despite having the 
same v r distribution. 

In the isolated thin disk (upper middle panel of Figure |5), the 
range of J z grows relative to the initial values at every energy E. 
The change is largest at lower energies. Recall from Figure [T]that 
the isolated thin disk develops a bar, which is composed of particles 
on radial orbits with relatively low J z , in the most bound region 
of the disk. Thus, this relatively large change towards less circu- 
lar orbits at low energies can be associated with bar formation. In 
contrast, the distribution of J z in the isolated thick disk is basically 
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Figure 4. The fraction of particles that migrate more than 1.0 kpc (solid 
line), 2.0 kpc (long dashed line), and 3.0 kpc (dotted line) as a function of 
formation radius. Results are binned such that the migration probability is 
calculated for particles in non-overlapping 1.0 kpc wide annuli. Lines con- 
nect the migration fraction in each bin (x-coordinates are bin centers, from 
5.5 to 19.5 kpc). The gray regions are the radial mass profiles of the ini- 
tial disks normalized such that the total mass contained in the first annulus 
equals 2/3 on this scale. The migration probability follows the mass distri- 
bution in the isolated disks but is anti-correlated with mass in the perturbed 
disks. 



equivalent to its initial state. In the isolated thick disk there is no 
bar, relatively little spiral structure develops, and we find the least 
radial migration among the four simulations. Ignoring the region of 
the Lindblad diagram influenced by members of the bar, the iso- 
lated disk diagrams suggest that the extensive radial migration seen 
in the thin disk and associated with guiding center modification de- 
creases the angular momentum of particles. 

The two perturbed disks (right column of Figure [5} show 
larger changes in their Lindblad diagrams. Bar formation in both 
simulations can explain the significantly lower median J z at lower 
energies. At higher energies, corresponding to less bound particles 
further out in the disk, both disks show a much larger dispersion 
in Jz than their isolated counterparts. Additionally, there is sub- 
structure in the Lindblad diagrams (groups of relatively low J z in 
a narrow range of energy) not seen in the isolated disks. Satellite 
bombardment in the perturbed disk simulations has a qualitatively 
discernible impact on the angular momentum distribution of the 
disk. 

To quanti fy changes in J z , we introduce the circularity (e) 
quantity (e.g., lAbadi et"al]|2003l) . For a particle of energy Ei and 
angular momentum Ji in the ^-direction, we define circularity as 
e — Ji / Jcirc (Ei), the ratio of Ji to the specific angular momentum 
the particle would have if it were on a circular orbit with energy Ei 
(obtained using the circular orbit curve). Circular orbits have e = 1 
and radial orbits have e = 0. Negative circularities correspond to 
retrograde o rbits. We no t e that this method is formally different 
from that of lAbadi et all J2003I) . The circular orbit curve may not 
represent the highest J z for a given E due to shot noise in the po- 
tential or gravitational potential asym metries. Strictly defi ning the 
maximum J z in each energy bin (as in lAbadi et al l d2003l) ) results 
in the same qualitative trends discussed later. Our method benefits 
from the use of a mathematically constructed and reproducible ro- 
tation curve and systematically decreases Ae= ti — ei at the 0.01 
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Figure 5. Lindblad diagrams for all collisionless simulations. Each column shows the initial (left); final, isolated (middle); and final, perturbed (right) particle 
angular momentum projected along the z axis (J z , kpc kms -1 ) vs. specific energy (E, km 2 s~ 2 ). The red, yellow, and orange regions encompass 68%, 
95%, and 99%, respectively, of particles centered on the median J z as a function of E. J z > 99% outliers are plotted as points. The green line denotes 
the theoretical angular momentum - energy curve for circular orbits. Each panel includes a histogram illustrating the distribution of circularity (e) for all disk 
particles (see text for details). The isolated thin disk shows evolution towards slightly non-circular orbits. There is little change in e in the isolated = 400 pc 
case. Simulations with satellite bombardment show a pronounced shift towards less circular orbits. For reference, the total energy of midplane orbits at 2, 5, 
fO, and f 5 kpc from the galactic center is approximately —2.0, —1.6, —1.2, and —1.0 km 2 s — 2 , respectively. 



level. The inset of each panel in Figure [5] shows the circularity dis- 
tribution of each simulation. 

Figure [6] plots the relation between circularity and eccentric- 
ity for orbits near the solar radius (8.0 kpc) in the Zd = 400 pc 
disk. Circularity is related to eccentricity via the particle's radial 
and vertical energies as well as its orbital inclination (since J z is 
a projected quantity). Orbits in the midplane of the disk have the 
highest circularity for a given eccentricity. Shaded regions in Fig- 
ure^ show the 95 th and 99* h percentile circularity as a function 
of eccentricity at the solar radius. For reference in interpreting Fig- 



ure [5] and subsequent figures, it is worth noting that changes of 
~ 0.02 in e typically correspond to quite noticeable changes in ec- 
centricity. The boundaries of these regions are not smooth due to 
binning and small number statistics for initially high eccentricity 
particles. 

Returning to Figure|5] we see that the two isolated disks show 
markedly different evolution of their circularity distributions. In 
the isolated thin disk, the fraction of stars with e ~ 1 (rightmost 
bin) drops from 0.38 to 0.26, and the median e drops from 0.980 
to 0.955. Particles in the bar predominantly populate the newly 



Radial Mixing in Galactic Disks 9 




0.9 - 



0.05 0.1 0.15 0.2 0.25 

eccentricity 

Figure 6. Possible values of circularity as a function of eccentricity in the 
solar annulus (7 ^ r ^ 9 kpc) of the jz<j = 400 pc disk. The top edge 
of the dark region represents orbits in the midplane of the disk. The dark 
region encompasses the 95 percentile of particle circularity in the solar 
annulus of each simulation; the bottom edge of the light region corresponds 
to the 99 th percentile circularity. 



formed low circularity tail. In the thick isolated disk, on the other 
hand, the e distribution is nearly identical to the initial disk's, with 
median e dropping only 0.002. This lack of evolution in the cir- 
cularity distribution is fully consistent with CR scattering. As de- 
scribed by SB02, individual stars may exchange angular momen- 
tum across CRs (churning) while the overall distribution would re- 
main unchanged. However, our results from Section UTI show that 
the Delta R distribution of this model, which reflects individual par- 
ticles and their orbits, is nearly identical to that expected from ob- 
serving the particles' initially elliptical orbits. The circularity dis- 
tribution of the isolated, thick disk combined with its Delta R dis- 
tribution imply that, on average, individual particle guiding centers 
are not significantly modified. 

The circularity distributions of the perturbed disks are demon- 
strably altered from their isolated counterparts but are similar to 
one another. The most common circularities are now in the range 
0.96 < e < 0.98 in both disks, with a decrease for e> 0.98. The 
thin disk starts with more e ~ 1 particles because of its smaller 
orbital inclinations, and its circularity distribution evolves more 
strongly, with median e dropping from 0.980 to 0.936 vs. 0.963 
to 0.929 for the perturbed thick disk. Notably, the perturbed thin 
and thick disks evolve to similar Ar (Figure |3) and e distributions 
despite starting with different scale heights. 

Figure [7J tracks the changes of selected individual particles in 
the (E, J z ) space of the Lindblad diagram. While E and J z are 
not individually conserved in the presence of a non-axisymmetric 
perturbation, the Jacobi invariant / = E — Qj, J z is, where Qt 
is the pattern speed of the pe rturbation, assumed to be static and 
small (SB02: ISellwoo3l2010h . If AI = 0, then AJ Z /AE « fl b . 
The SB02 mechanism operates at the corotation resonance of the 
star/particle and the spiral wave, requiring that fl;, = f2 ro t- Thus, 
in the galaxy, particles should move parallel to the line that is tan- 
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Figure 7. The initial and final E, J z pairs for ten randomly selected parti- 
cles in each simulation with initial energy —1.405 E[ ^ —1.395 X 10 
km 2 s — 2 and angular momentum 1.645 ^ J z ,i ^ 1.655 X 10 3 kpc 
kms -1 (large, filled square). For clarity we require that plotted particles 
lose at least 0.05 X 10 5 km 2 s -2 in energy during the simulation. The 
final E,J Z pairs of the ten particles are plotted as open squares. Lines con- 
nect the initial and final E, Jz points of each particle. The E, J z curve 
populated by circular orbits in the initial state of each simulation is indi- 
cated by the thick black curve. The dashed line is tangent to the circular 
orbit curve at the initial energy of all ten particles. 



gent to the circular orbit curve at their binding energy prior to scat- 
tering. In other words, the SB02 mechanism requires that changes 
in energy be accompanied by changes in angular momentum that 
preserve the orbital shape, modulo differences in the slope of the 
circular orbit curve over the range [Ei, Ef] of a given particle. 

Figure [7] zooms in on the area of the Lindblad diagram desig- 
nated by the two small boxes drawn on the initial state diagrams in 
Figure[5] We randomly select ten particles within the same initial E 
and J z range (filled, dark squares in Figure [7} and plot their final E 
and J z . We require Et-Ei = AE < -0.05xl0 5 km 2 s~ 2 toen- 
sure that the individual tracks are visible. In the isolated thin disk 
(top left), most particles move parallel to the circular orbit curve 
tangent (dashed line); this relationship between AJ Z and AE is 
consistent with corotational resonant scattering (equation 2, SB02). 
When particles lose energy in the isolated thick disk (bottom left), 
their J z typically remains closer to the circular orbit curve than 
in the isolated thin disk. The changes in energy and angular mo- 
mentum are relatively small, indicating that the guiding centers are 
modified to a modest extent dBinnev & Tremaine]|2008h . Note that 
these particles are in the top 1% of |AJ Z |; most particles in this 
particular experiment have A J- ~ 0. The randomly selected parti- 
cles in both perturbed disks have slopes AJ Z / AE that are steeper 
than the slope of the tangent to the circular orbit curve at corota- 
tion. This distinct coupling of AE and AJ Z , combined with our 
results concerning the migration probability as a function of Rf ornl 
(Section 13. It . offer compelling evidence that migration in the per- 
turbed disks is driven, at least in part, by a mechanism that does not 
operate in the isolated systems. Figure[7Jshows that the orbital char- 
acteristics of many particles in the perturbed disks are modified in 
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Figure 8. Median change in circularity Ae as a function of Ar for the entire 
disk (open squares) and for particles with rf , n > 4.0 kpc (filled squares). 
Particles are sorted into non-overlapping 0.5 kpc bins of Ar. Error bars 
mark the 25 th and 75 th percentile Ae in each bin. Histograms at the bot- 
tom of each panel represent the relative fraction of particles in each bin for 
the case rf , r; > 4.0 kpc. 



a fashion inconsistent with either epicyclic motion or a single spiral 
wave scattering event. 

We now examine the correlation between migration and 
changes in orbital properties. Since the metallicity of star-forming 
gas increases with time and decreases with radius, any such corre- 
lations also imply observable correlations between the present or- 
bital parameters and metallicities of stars as a function of age and 
Galacto-centric radius. Figure[8]shows the median change in circu- 
larity, Ae = ef — ei, as a function of radial migration distance Ar 
in the four collisionless simulations (open squares). In the isolated 
thick disk; changes in circularity are tiny (median |Ae| ^ 0.005) 
at any Ar; Figure [7J shows that particles in this simulation stay 
close to the circular velocity curve even if they change energy. In 
the other three simulations, particles with negative Ar show sub- 
stantial drops in circularity (typical median Ae < —0.06), which 
are almost certainly associated with bar formation. The bars that 
form in these three simulations increase the orbital eccentricities of 
their members. To focus on behavior in the disk proper, the filled 
squares in each panel show the median Ae vs. Ar for those par- 
ticles that start and end the simulation at r ^ 4 kpc, beyond the 
extent of the bar. Error bars mark the inter-quartile range (25 th to 
75 th percentile) at each Ar. 

In the isolated thin disk, particles that migrate inwards (Ar< 
0) experience a modest decrease in circularity, stronger for more 
negative Ar (and AE), consistent with the tracks shown in Fig- 
ure [JJ Particles that migrate outwards (Ar> 0) increase their en- 
ergy and typically traverse areas of the Lindblad diagram with rel- 
atively little change in the slope of the circular orbit curve. Fol- 
lowing the tangent to the circular orbit curve, particles will not 
significantly change their circularity in such a scenario (note the 
relative lack of low circularity particles at E > —1.5 km 2 s -2 in 
Figure [5). The perturbed disks show a decrease in median circular- 
ity at every Ar, and a much wider inter-quartile range indicating 
a greater range of orbital inclinations and eccentricities . Circular- 
ity drops more strongly for particles that have experienced strong 
radial migration, either inward or outward. The minimum in these 
curves is slightly offset to positive Ar because dynamical heat- 
ing slightly puffs up the disk radially, decreasing the potential at 
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Figure 9. Like Figure[8] but the change in maximum vertical displacement, 
Az max , is plotted in place of the change in circularity. Open squares show 
all disk particles, filled squares show those with rf , r; > 4 kpc, and error 
bars mark 25 th and 75 th percentile at a given Ar. 



a given radius, thereby increasing its total energy, moving parti- 
cles to the right on the Lindblad diagram , and lowering the circu- 
larity for particles with Ar= 0. Figure [8] indicates that stars with 
anomalous chemistry for their age and current position should have 
preferentially more eccentric orbits. The trend is smaller than the 

inter -quartile range, but similar in mag nitude. 

ISchonrich & Binnevl d2009h and lLoebman et"al] | |2010|) dis- 
cuss the possible role of radial migration in producing thick d isks 
like the ones observed in the Milky Way dGilmore et al.lll989l) and 
other edge-on galaxies jDalcanton & Bernstein! 20021) . Figure [9l is 
similar to Figure [8] but instead of Ae it plots the change in verti- 
cal energy, quantified by the maximum distance z max that a par- 
ticle can reach from the disk plane. We compute z max approx- 
imately from each particle's vertical velocity component assum- 
ing that the final potential of each simulation is static, namely 



4jrGE(r) 



where z is the vertical position of the parti- 



cle at the final output, v z is the velocity along the z axis, G is the 
gravitational constant, and E(r) is the surface density of the disk 
at radius r assuming all the mass of the disk is in the midplane. 
The change in z max is Az max = z max ,f — z maXii where the sub- 
scripts i and / refer to the initial and final simulation snapshots, 
respectively. 

The two isolated simulations show a shallow linear trend be- 
tween Az max and Ar (open squares for all particles; filled squares 
for those with n,r(> 4 kpc). When particles move outwards (in- 
wards) through the disk, they experience a weaker (stronger) grav- 
itational potential, thus increasing (decreasing) z max . However, 
changes in z max are small, less than 300 pc even when we con- 
sider the quartile range at the extremes of Ar. The perturbed disks, 
by contrast, show larger changes in the median z max and a dramatic 
increase in the inter-quartile range at all Ar. K08 show that the per- 
turbed 400 pc disk develops a two-component vertical structure in 
quantitative agreement with the observed thin/thick disk structure 
of the MW. The perturbed 200 pc disk increases its scale height 
(to ~ 500 pc at the solar annulus), but it can still be described by 
a single component model. This vertical heating by satellite per- 
turbations is evident in Figure [9] Particles that have large positive 
Ar have the largest increase in z max , which is plausibly a conse- 
quence of moving outwards to regions of lower disk surface density 
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Figure 10. The distribution of formation radius (-Rf or m. m kpc) for all 
particles that are in the solar annulus (7 kpc ^ r ^ 9 kpc) at the end of 
each simulation. Histograms report the fraction of solar annulus particles 
emigrating from non-overlapping 500 pc annuli in i?f or m ■ Particles within 
the dashed lines remained in the solar annulus throughout the simulation. 
Both thin disks and the perturbed thick disk show a broad range of iJf orm 
at the solar annulus. The gray histogram in each panel is the Ri OTm distri- 
bution for those particles that are in the solar annulus and at large heights 
above the plane (1.0 < \z\ < 1.5 kpc) in the final simulation output. 

and thus weaker vertical restoring force. For Ar < 0, the trend of 
median Az max with Ar is approximately flat, suggestion a cancel- 
lation between the effects of increased E(r) at smaller r and direct 
satellite-induced heating of those particles with the largest excur- 
sions. Figure [9]implies that stars with high metallicity for their age 
and present location should have preferentially larger z max , though 
the scatter is larger than the trend. 

3.3 Solar Annulus 

The solar neighborhood is easier to study than other regions of 
the Galaxy, since high-precision spectroscopy is easier for brighter 
stars and parallax and proper motion measurements are more ac- 
curate at smaller distances. Some of the most detailed chemo- 
dynamic surveys, such as the GCS dNordstrom et al .1120041) and the 
Radial Velocity Experiment (RAVE, ISteinmetal2006l) concentrate 
on the solar neighborhood. In this section, we repeat some of our 
earlier analysis specifically for stars that reside in the solar annulus 
(7 kpc ^ rf 9 kpc) at the end of each simulation. This focus on 
the solar annulus also removes much of the impact of the bars that 
dominate evolution of the inner disk (r < 3 kpc) in three of our 
simulations, though some particles from the bar region can migrate 
as far as the solar radius, and resonances betwe en the bar and spiral 
struct ure may increase migration frequency dMinchev & Famaevl 
l201fjh . 

Figure [TO] shows the radius of formation (i?f or m) distribution 
of particles residing in the solar annulus, marked by the vertical 
dashed lines. Only the isolated thick disk simulation predicts a fi- 
nal solar annulus dominated by stars bom in the solar annulus, with 
tails extending 1-2 kpc on either side. The broad Ar distributions 
of the other three simulations show that their stars migrate to the 
solar annulus from a wide range of formation radii. The global 
Ar (Figure [3} and solar annulus 7if orm distributions of the iso- 
lated and perturbed thin disks are remarkably similar, despite the 
differences in migration mechanisms discussed in Section [3T2l In 
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Figure 11. The change in circularity Ae as a function of Ar for particles 
ending the four collisionless simulations in the solar annulus (7 kpc ^ 
rf < 9 kpc). Particles are sorted into non-overlapping 0.5 kpc bins of 
Ar. We plot the median (squares) and the 25 th and 75 th percentile (error 
bars) Ae in each bin. Histograms at the bottom of each panel represent the 
relative fraction of particles in each bin. 



the isolated thin disk, 32% of solar annulus particles originated at 
rJf orm ^ 6 kpc and 4% at i?f or m ^ 10 kpc. Corresponding numbers 
for the perturbed thin disk are 36% and 3%. The Rf OI m distribution 
of the perturbed thick disk is slightly narrower, but it still broad 
with respect to the isolated thick disk. 

Figure[TT]shows the correlations between Ae and Ar for par- 
ticles that end in the solar annulus. Consistent with results for the 
full disk (Figure[8j, the solar annulus particles in the isolated thick 
disk show no significant change in median e regardless of Ar. In 
the isolated thin disk, the range of Ae is much larger, with a mod- 
est decrease in median e. The median Ae drops at large positive Ar 
because particles move to less inclined orbits as they migrate out- 
wards and disk heating has slightly modified the galaxy's circular 
velocity curve, allowing outward moving particles to potentially in- 
crease circularity. In the perturbed thick disk, there is a strong and 
nearly linear trend between Ae and Ar. Particles that migrated to 
the solar annulus from the outer disk have experienced substantial 
drops in circularity, while particles migrating from the inner disk 
show only modest decreases. The range of Ae is large at all Ar. 
Results for the perturbed thin disk are intermediate between those 
of the isolated thin and perturbed thick disks: quasi-linear trends at 
large | Ar j but a flat plateau at intermediate j Ar]. This behavior is 
plausibly a consequence of two different mechanisms contributing 
to migration, with spiral-induced mixing dominating at intermedi- 
ate Ar and satellite-induced mixing dominating at the extremes. 

Figure[J_2]shows the solar annulus correlations of Az max with 
Ar, analogous to Figure [8] for the full disk. For the two isolated 
disks there is a clear linear trend of median Az max with Ar as 
expected from the arguments in Section IT21 particles that migrate 
outward move to a region of lower disk surface density, so if their 
vertical velocities are not systematically changed by the radial mi- 
gration they will attain higher 2 max . Both perturbed disks show 
signs of the strong vertical heating induced by satellite accretion. 
The median Az max is ^ 0.2 kpc almost independent of Ar, and 
the range of Az max is much larger than in the isolated disks. Since 
inwardly migrating particles experience a higher vertical potential 
at rf than n, they must experience more vertical heating than out- 
wardly migrating particles to keep the Az max trend flat. 
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Ar [kpc] 

Figure 12. Like Figure [TTI but the change in vertical displacement Az n 
is plotted in place of the change in circularity. 



Returning to Figure [TU] gray histograms represent the i?f or m 
distributions of particles that end the simulations in the solar an- 
nulus at high z, 1.0 < z < 1.5 kpc. In the isolated thick disk 
and both perturbed disks, the -Rform distribution of high z particles 
resembles that of all solar annulus particles. In these three models, 
therefore, selecting high- z particles does not isolate a population 
with atypical radial migration. In the isolated thin disk, on the other 
hand, the i?f or m distribution of high-z particles is strongly skewed 
towards low formation radii, with a peak at -Rf rm= 5 kpc. Here, 
only initially hot particles that migrate significantly outwards via 
churning or interactions with the bar and experience a weaker po- 
tential have enough vertical energy to overcome the local restoring 
force and reach high z. Only 0.1% of the solar annulus is at high 
z in the isolated thin disk (truly the tail of the initial vertical ve- 
locity dispersion); this number rises to 3.2% and 7.5% in the thin 
and thick perturbed disks, respectively. The radial migration mech- 
anisms in the perturbed disks ensure that even the high z population 
of the solar annulus comes from a broad range of Rform- Measure- 
ments of the age-metallicity relation for high-z stars could be a 
valuable diagnostic for distinguishing models of radial migration 
and vertical heating. 



4 SUMMARY AND DISCUSSION 

Observations suggest that radial migration pla ys an important role 
in the chemical evolution of the MW disk dWielenetal1ll996l : 
ISchonrich & BinnevkOO !) . SB02 described a mechanism by which 
stars at corotational resonance with spiral waves can scatter, pro- 
ducing large changes in guiding center radius while keeping stars 
on nearly circular orbits. Previous simulations have shown that mi- 
gration over several kiloparsecs can occu r in isolated disks grown 
by smooth accretion dRoskar et al.l2008a | ) and that encou nters with 
satellites can also induce migration dOuillen et al1l2009l) . Here we 
have carried out a systematic investigation of a variety of simula- 
tions to characterize the role of stellar disk properties, gas fractions, 
and satellite perturbations in producing radial migration. Most im- 
portantly, our suite of simulations includes experiments with a level 
of satellite bombardment expected in ACD M models of galaxy for - 
mation, which earlier investigations (K08, iKazantzidis et al.l2009h 
have shown to produce vertical and in-plane structure resembling 
that seen in the MW. 



For disks evolved in isolation, the degree of migration cor- 
relates with the degree of disk self-gravity, susceptibility to bar 
formation, and spiral structure. The collisionless stellar disk with 
Zd = 400 pc, chosen to match that of the thin disk in the MW, 
has no bar and minimal spiral structure. It exhibits limited radial 
migration; the median value of radial change \n — n is Ar mc( j = 
0.47 kpc, consistent with the level expected from epicyclic mo- 
tion (Figure [3}- The collisionless Zd = 200 pc disk is much more 
unstable, develops a bar and spiral structure (Figure [T], and has 
a higher degree of radial migration, with Armed = 0.91 kpc and 
Ar8o= 2.03 kpc. The presence of gas is a catalyst for radial mi- 
gration in the Zd = 400 pc case. Both Ar mc d and Ar$o increase in 
models with larger gas fractions or lower star formation efficiency 
(which consumes gas more slowly). The galactic bar strongly in- 
fluences, and perhaps dominates, the resultant individual particle 
dynamics in the inner galaxy. Substantial radial migration, con- 
sistent with the SB02 mechanism or being satellite-induced, oc- 
curs in the outer portion of the disk in all of our collisionless 
experiments except for the isolated Zd = 400 pc disk. In the 
Zd = 200 pc disk, stars that finish the simulation in the solar annu- 
lus (7 < R < 9 kpc) have Ar mcd = 1.36 kpc, Ar 8 o= 2.68 kpc. 

Adding satellite perturbations dramatically increases radial 
migration in the Zd = 400 pc disks, raising Ar mc d from 0.47 kpc 
to 0.89 kpc globally and from 0.60 kpc to 1.15 kpc in the solar 
annulus. Perturbations do not change the Ar distribution so drasti- 
cally in the Zd = 200 pc disks, but there are other indications that 
the nature of radial migration is different. In the isolated galaxies, 
migrating particles in the outer disk follow tracks in E, J z space 
that parallel the tangent to the local circular velocity curve, as ex- 
pected for the SB02 mechanism (see Section However, in- 
wardly migrating particles in the perturbed disks lose more angular 
momentum for a given change in energy, causing migration asso- 
ciated with heating and thus different from the SB02 mechanism 
in this respect. The relationship between AE and AJ Z found in 
the perturbed experiments is broadly more consistent with scatter- 
ing at LR (see SB02), but a detailed decomposition of the modes 
in the disk and satellites (beyond the scope of this paper) is neces- 
sary to confirm LR scattering on a particle by particle basis. The 
radial distribution of migrating particles presents an even clearer 
distinction, one with important implications for chemical evolu- 
tion models (Figure [4). In the isolated disks, the probability that 
a particle with a formation radius -Rf or m undergoes significant mi- 
gration ( | Ar > 1 kpc) is approximately proportional to the disk 
surface density at Rt om ; i.e., migration follows mass. In the per- 
turbed disks, the probability of migration is flat or increasing with 
radius, so a much larger fraction of migration comes from the tenu- 
ous outer disk. Particles in the outer disk are overall more likely to 
migrate relative to the inner disk as 1) they feel a weaker potential 
and are thus more susceptible to direct heating by satellites and 2) 
their circular frequencies make them more likely to resonantly in- 
teract with satellites and migrate. Thus, satellite-induced migration 
patterns are distinct from those produced by purely secular evolu- 
tion. 

Satellite bombardment heats the disk both vertically and ra- 
dially. Compared to the isolated disks, the perturbed disks experi- 
ence a greater drop in median circularity (hence growth of median 
eccentricity) during their evolution, and individual particles experi- 
ence a wider range of circularity changes. The circularity change is 
moderately correlated with radial migration distance and direction, 
though the trend is smaller than the inter-quartile range at a given 
Ar (Figure |8). If we restrict our analysis to the solar annulus, the 
correlation between Ar and Ae is more significant. In particular, 
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particles in the solar annulus of the = 400 pc disk that migrate 
2-4 kpc inwards experience substantial drops in circularity, while 
those that migrate the same distance outwards nearly maintain their 
initial circularity (Figure [TTV The large inward excursions from the 
outer disk, driven by satellite perturbations, systematically remove 
angular momentum from particle orbits. 

As shown by K08, vertical heating by satellite bombardment 
produces (in the case of the z& = 400 pc disk) a two-component 
vertical structure in quantitative agreement with the thin and thick 
disk profiles observed in the MW. In the isolated disks, we find 
the expected trend that as particles migrate outwards, they expe- 
rience a weaker vertical potential and increase their vertical en- 
ergy (characterized by 2 max , the maximum distance from the plane 
that a particle can reach given its current location and velocity). 
However, in our simulations, this effect is not sufficient to produce 
a second, "thick-disk" component in the isolated systems even if 
they have substantial radial migration. The satellite perturbations 
increase the median and range of Zmax considerably at all Ar. The 
resulting overall trend of Az max with Ar is fairly flat (especially at 
a fixed r{ , such as the solar annulus, see Figuresl9landll2t. suggest- 
ing that systems with low vertical velocity dispersion and subjected 
to weaker potentials (as in the outer disk) are more susceptible to 
vertical heating (similar to the radial heating trends seen in Sec- 
tion l3.lt . Particles with -Rf or m interior or exterior to a given rf can 
be found at significant distances above the plane in the perturbed 
disks, while increases in 2 max are smaller in the isolated systems 
and rely on particles moving outwards (Figure HO], 

Our results confirm earl ier findings that spiral structure de- 
velop ment in isolated disks l lRoskar et al.ll2008al. iLoebman et all 



|2010|) and perturbations by satellites dOuillen et al.ll2009f) can pro- 
duce significant radial mixing of stellar populations while retaining 
reaso nable orbital structu re for disk stars. They strongly sup port the 
view Jwielenet all 1994 SB02: ISchonrich & Binnevl2009n that ra- 
dial mixing is an essential ingredient in understanding the chemical 
evolution of the MW and disk galaxies in general. While a combi- 
nation of metallicity and age can be used to estimate a star's forma- 
tion radius given a chemical evolution model, the uncertainties in 
this approach (including the difficulty of estimating ages for typical 
stars in a spectroscopic survey) will make it difficult to reconstruct 
formation radius distributions for observed stars, even in the solar 
neighborhood. However, our analysis provides theoretical guidance 
for chemical evolution mo dels that incorporate radial mixing (e.g. 
ISchonrich & Binnevl2009l) and suggestions for the correlations one 
might search for between chemical abundances and orbital proper- 
ties (though the predicted trends are fairly weak). We regard our 
perturbed z& = 400 pc disk as the most relevant simulation, as it 
includes the satellite accretion events expected in ACDM and pro- 
duces a final vertical structure similar to that measured in the MW. 
In this simulation 41% of particles in the 7 kpc < R < 9 kpc so- 
lar annulus were "born" in that annulus, 20% migrated there from 
R < 6 kpc, and 7% migrated there from R > 10 kpc. 

Radial mixing and orbital dynamics changes are sensitive to 
several different aspects of disk modeling, as shown by the vari- 
ety of our results. Robust predictions should therefore come from 
calculations that self-consistently include star formation, chemical 
enrichment, gas accretion, and accretion events, and the interplay 
among these elements. We will move in this direction with our 
future simulations. Giant spectroscopic surveys such as SEGUE, 
RAVE, APOGEE, LAMOST, and HERMES offer an extraordinary 
opportunity to unravel the formation history of the MW, and they 
offer an exciting new challenge to theoretical models of galaxy for- 
mation. 
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